Metabarcoding Singapore
ASV from Dada2
Metabarcoding Singapore
ASV from Dada2
- 1 Aim
- 2 Analysis
- 2.1 Assignment of eukaryotic ASVs based on PR2 database
- 2.2 Phyloseq analysis
- 2.2.1 Create phyloseq files for euk after filtering the data
- 2.2.2 Create phyloseq files for proks
- 2.2.3 Create a list for the auto and hetero phyloseq files
- 2.2.4 Bar plot of divisions per station for auto and hetero
- 2.2.5 Bar plot of class per station for auto and hetero
- 2.2.6 Compare by Straight, Site, Moonsoon (abundant OTUs only)
- 2.2.7 Main genera for different division of Eukaryotes (Autotrophs)
- 2.2.8 Main species of Mamiellophyceae and Diatoms (Eukaryotes - Autrotrophs)
- 2.2.9 Heatmaps (abundant OTUs)
- 2.2.10 NMDS
- 2.2.11 Network analysis
Initialize. This file defines all the necessary libraries and variables
source('Metabarcoding Singapore_init.R', echo=FALSE)1 Aim
- Assign and analyze eukaryotes for Singapore metabarcoding data (ASV assigned with dada2 as implemented on Mothur).
- Do some analyzes with the prokaryotes too…
2 Analysis
2.1 Assignment of eukaryotic ASVs based on PR2 database
2.1.1 File names
full_path_data <- function(file_name) {
str_c("../qiime/2018-09-06_dada2/", file_name)
}
taxo_file <- full_path_data("taxonomy.tsv")
otu_file <- full_path_data("feature-table_unrarefied.tsv")
sequence_file <- full_path_data("ref_sequences_unrarefied.fasta")
metadata_file <- "../metadata/Singapore_metadata.xlsx"
sequence_file_euk <- full_path_data("ASV_unrarefied_euk.fasta")
dada2_taxo_file_euk <- full_path_data("ASV_unrarefied_euk.dada2.taxo")
dada2_boot_file_euk <- full_path_data("ASV_unrarefied_euk.dada2.boot")
otu_table_final_file <- full_path_data("ASV_final.tsv")
blast_file <- full_path_data("ASV_unrarefied_euk.blast.tsv")2.1.2 Read the data
- The dada2 treatment has already removed the forward and reverse primers, so no need to remove them
- Work with the unrarefied data
# Read the sample and metadata tables
sample_table <- read_excel(metadata_file, sheet="samples", range="A1:C89")
metadata_table <- read_excel(metadata_file, sheet="metadata", na=c("ND", ""))
sample_table <- left_join(sample_table, metadata_table) %>%
mutate(sample_label = str_c(strait_label,location_label,
monsoon,sprintf("%03d",day_number),
sep="_"))
# Read the taxonomy table
taxo_table <- read_tsv(taxo_file)
# Clean up the taxonomy
taxo_table <- taxo_table %>%
mutate(taxo_clean = str_replace_all(Taxon, "D_[0-9]+__","")) %>%
separate(col=taxo_clean, into=str_c("taxo", c(1:7)), sep=";") %>%
rename(otu_name = `Feature ID`)
# Read the otu table
otu_table <- read_tsv(otu_file, skip=1) %>% # Jump the first line
rename(otu_name = `#OTU ID`) %>%
mutate(otu_id = str_c("otu_", sprintf("%04d",row_number())))
# Read the sequences
otu_sequences <- readAAStringSet(sequence_file)
otu_sequences.df <- data.frame (otu_name=names(otu_sequences),sequence=as.character(otu_sequences))
# Remove the primers - Not necessary because the primers have been removed
# fwd_length = 20
# rev_length = 15
# otu_sequences.df <- otu_sequences.df %>%
# separate (col=names, into=c("otu_id_qiime", "otu_rep_seq"), sep=" ") %>%
# mutate (sequence = str_sub(sequence, start=fwd_length+1, end = - rev_length - 1))
otu_table <- taxo_table %>%
left_join(otu_table) %>%
left_join(otu_sequences.df) %>%
arrange(otu_id)
# Write a fasta file for blast with all taxonomy roups
# otu_sequences <- otu_table %>% transmute(sequence=sequence, seq_name=otu_id)
# fasta_write(otu_sequences, file_name="../qiime/otu_rep_98_all.fasta", compress = FALSE, taxo_include = FALSE) 2.1.3 Only keep the eukaryotes in the OTU file
otu_table_euk <- otu_table %>% filter(str_detect(Taxon, "Eukaryota"))
# Write the fasta file file
otu_sequences_euk <- otu_table_euk %>% transmute(sequence = sequence, seq_name = otu_id)
fasta_write(otu_sequences_euk, file_name = sequence_file_euk, compress = FALSE,
taxo_include = FALSE)[1] TRUE
2.1.4 Use dada2 to reassign to PR2
dada2_assign(seq_file_name = sequence_file_euk)2.1.5 Read the PR2 assignement and merge with initial otu table
otu_euk_pr2 <- read_tsv(dada2_taxo_file_euk)
otu_euk_pr2_boot <- read_tsv(dada2_boot_file_euk) %>% rename_all(funs(str_c(.,
"_boot"))) %>% rename(seq_name = seq_name_boot)
otu_euk_pr2 <- left_join(otu_euk_pr2, otu_euk_pr2_boot) %>% rename(otu_id = seq_name)
otu_table_final <- left_join(otu_table, otu_euk_pr2) %>% select(otu_id, otu_name,
taxo1:taxo7, Taxon, kingdom:species_boot, matches("EC|PR|RM|SBW|STJ"), sequence)
write_tsv(otu_table_final, otu_table_final_file, na = "")2.1.6 Process BLAST file
blast_18S_reformat(blast_file)2.2 Phyloseq analysis
2.2.1 Create phyloseq files for euk after filtering the data
Filter the euk data to remove the low bootstraps values (threshold : bootstrap > 90% at the supergroup level) and create a phyloseq file
Note the bootstrap threshold had to be higher for 98% compared to 97% (90% vs 65%)
otu_table_euk_final <- otu_table_final %>% filter(supergroup_boot > 90)
str_c("Number of euk otus : ", nrow(otu_table_euk_final))[1] "Number of euk otus : 919"
otu_mat <- otu_table_euk_final %>% select(otu = otu_id, matches("EC|PR|RM|SBW|STJ"),
-species, -species_boot)
tax_mat <- otu_table_euk_final %>% select(otu = otu_id, kingdom:species)
samples_df <- sample_table %>% rename(sample = sample_id)
row.names(otu_mat) <- otu_mat$otu
otu_mat <- otu_mat %>% select(-otu)
row.names(tax_mat) <- tax_mat$otu
tax_mat <- tax_mat %>% select(-otu)
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select(-sample)
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
singa_euk <- phyloseq(OTU, TAX, samples)
singa_eukphyloseq-class experiment-level object
otu_table() OTU Table: [ 919 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 919 taxa by 8 taxonomic ranks ]
Break up into photosynthetic and non-photosynthetic (Metazoa are removed)
singa_photo <- subset_taxa(singa_euk, (division %in% c("Chlorophyta", "Dinoflagellata",
"Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta", "Cercozoa")) &
!(class %in% c("Syndiniales", "Sarcomonadea")))
singa_photophyloseq-class experiment-level object
otu_table() OTU Table: [ 312 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 312 taxa by 8 taxonomic ranks ]
singa_hetero <- subset_taxa(singa_euk, !(division %in% c("Chlorophyta", "Dinoflagellata",
"Cryptophyta", "Rhodophyta", "Haptophyta", "Ochrophyta", "Cercozoa", "Metazoa")) |
(class %in% c("Syndiniales", "Sarcomonadea")))
singa_heterophyloseq-class experiment-level object
otu_table() OTU Table: [ 366 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 366 taxa by 8 taxonomic ranks ]
Normalize number of reads in each sample using median sequencing depth.
total_photo = median(sample_sums(singa_photo))
sprintf("The number of reads used for normalization of autotrophs is %.0f",
total_photo)[1] "The number of reads used for normalization of autotrophs is 3253"
standf = function(x, t = total_photo) (if (sum(x) > 0) {
t * (x/sum(x))
} else {
x
})
singa_photo = transform_sample_counts(singa_photo, standf)
total_hetero = median(sample_sums(singa_hetero))
sprintf("The number of reads used for normalization of heterotrophs is %.0f",
total_hetero)[1] "The number of reads used for normalization of heterotrophs is 777"
# ! If there no cells do not transform
standf <- function(x, t = total_hetero) (if (sum(x) > 0) {
t * (x/sum(x))
} else {
x
})
singa_hetero = transform_sample_counts(singa_hetero, standf)Remove taxa which are in low abundance (< 0.10 in any given sample) and normalize again…
contrib_min <- 0.1
# Photo abundant
singa_photo_abundant <- filter_taxa(singa_photo, function(x) sum(x > total_photo *
contrib_min) > 0, TRUE)
total = median(sample_sums(singa_photo_abundant))
sprintf("The number of reads used for normalization of abundant autotrophs is %.0f",
total)[1] "The number of reads used for normalization of abundant autotrophs is 2946"
standf = function(x, t = total) (t * (x/sum(x)))
singa_photo_abundant = transform_sample_counts(singa_photo_abundant, standf)
singa_photo_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 65 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 65 taxa by 8 taxonomic ranks ]
# Hetero abundant
singa_hetero_abundant <- filter_taxa(singa_hetero, function(x) sum(x > total_hetero *
contrib_min) > 0, TRUE)
total = median(sample_sums(singa_hetero_abundant))
sprintf("The number of reads used for normalization of abundant heterotrophs is %.0f",
total)[1] "The number of reads used for normalization of abundant heterotrophs is 583"
standf = function(x, t = total) (t * (x/sum(x)))
singa_hetero_abundant = transform_sample_counts(singa_hetero_abundant, standf)
singa_hetero_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 91 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 91 taxa by 8 taxonomic ranks ]
2.2.2 Create phyloseq files for proks
otu_table_prok <- otu_table %>% filter(taxo1 %in% c("Bacteria", "Archaea"))
otu_mat <- otu_table_prok %>% select(otu = otu_id, matches("EC|PR|RM|SBW|STJ"))
tax_mat <- otu_table_prok %>% select(otu = otu_id, taxo1:taxo7) %>% rename(kingdom = taxo1,
supergroup = taxo2, division = taxo3, class = taxo4, order = taxo5, family = taxo6,
genus = taxo7) %>% mutate(species = NA)
samples_df <- sample_table %>% rename(sample = sample_id)
row.names(otu_mat) <- otu_mat$otu
otu_mat <- otu_mat %>% select(-otu)
row.names(tax_mat) <- tax_mat$otu
tax_mat <- tax_mat %>% select(-otu)
row.names(samples_df) <- samples_df$sample
samples_df <- samples_df %>% select(-sample)
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)
OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)
singa_prok <- phyloseq(OTU, TAX, samples)Normalize number of reads in each sample using median sequencing depth.
total_prok = median(sample_sums(singa_prok))
sprintf("The number of reads used for normalization of prokaryotes is %.0f",
total_prok)[1] "The number of reads used for normalization of prokaryotes is 54273"
standf = function(x, t = total_prok) (t * (x/sum(x)))
singa_prok = transform_sample_counts(singa_prok, standf)
singa_prokphyloseq-class experiment-level object
otu_table() OTU Table: [ 2079 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 2079 taxa by 8 taxonomic ranks ]
Remove otus which are in low abundace (< 0.05 in any given sample) and normalize again
contrib_min <- 0.05
singa_prok_abundant <- filter_taxa(singa_prok, function(x) sum(x > total_prok *
contrib_min) > 0, TRUE)
total = median(sample_sums(singa_prok_abundant))
sprintf("The number of reads used for normalization of abundant prokaryotes is %.0f",
total)[1] "The number of reads used for normalization of abundant prokaryotes is 24353"
standf = function(x, t = total) (t * (x/sum(x)))
singa_prok_abundant = transform_sample_counts(singa_prok_abundant, standf)
singa_prok_abundantphyloseq-class experiment-level object
otu_table() OTU Table: [ 44 taxa and 81 samples ]
sample_data() Sample Data: [ 81 samples by 21 sample variables ]
tax_table() Taxonomy Table: [ 44 taxa by 8 taxonomic ranks ]
2.2.3 Create a list for the auto and hetero phyloseq files
ps_list <- list(ps = c(singa_photo, singa_hetero, singa_prok), title = c("Eukaryotes - Autotrophs - all OTUs",
"Eukaryotes - Heterotrophs - all OTUs", "Prokaryotes - all OTUs"))
ps_list_abundant <- list(ps = c(singa_photo_abundant, singa_hetero_abundant,
singa_prok_abundant), title = c("Eukaryotes - Autotrophs - abundant OTUs (> 10%)",
"Eukaryotes - Heterotrophs - abundant OTUs (> 10%)", "Prokaryotes - abundant OTUs (> 5%)"))2.2.4 Bar plot of divisions per station for auto and hetero
Note: some stations are completely missing heterotrops (Only Opistokonta)
for (i in 1:3) {
p <- plot_bar(ps_list$ps[[i]], x = "sample_label", fill = "division") +
geom_bar(aes(color = division, fill = division), stat = "identity",
position = "stack") + ggtitle(str_c("Division level - ", ps_list$title[[i]])) +
theme(axis.text.y = element_text(size = 10)) + theme(axis.text.x = element_text(angle = 0,
hjust = 0.5)) + coord_flip()
print(p)
}2.2.5 Bar plot of class per station for auto and hetero
Only consider the abundant taxa
for (i in 1:3) {
p <- plot_bar(ps_list_abundant$ps[[i]], x = "sample_label", fill = "class") +
geom_bar(aes(color = class, fill = class), stat = "identity", position = "stack") +
ggtitle(str_c("Class level - ", ps_list_abundant$title[[i]])) + theme(axis.text.y = element_text(size = 10)) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5)) + coord_flip()
print(p)
}2.2.6 Compare by Straight, Site, Moonsoon (abundant OTUs only)
for (factor in c("strait", "location", "monsoon")) {
for (i in 1:3) {
ps_aggregate <- merge_samples(ps_list_abundant$ps[[i]], factor)
p <- plot_bar(ps_aggregate, fill = "division") + geom_bar(aes(color = division,
fill = division), stat = "identity", position = "stack") + ggtitle(str_c(ps_list_abundant$title[[i]],
" - ", factor)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5,
hjust = 1))
print(p)
}
}2.2.7 Main genera for different division of Eukaryotes (Autotrophs)
for (one_division in c("Chlorophyta", "Dinoflagellata", "Ochrophyta")) {
ps_subset <- subset_taxa(singa_photo_abundant, division %in% one_division)
p <- plot_bar(ps_subset, x = "genus", fill = "genus", facet_grid = strait ~
monsoon) + geom_bar(aes(color = genus, fill = genus), stat = "identity",
position = "stack") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5,
hjust = 1)) + ggtitle(str_c(one_division, " - Abundant OTUs"))
print(p)
}2.2.8 Main species of Mamiellophyceae and Diatoms (Eukaryotes - Autrotrophs)
for (one_class in c("Mamiellophyceae", "Bacillariophyta")) {
ps_subset <- subset_taxa(singa_photo_abundant, class %in% one_class)
p <- plot_bar(ps_subset, x = "species", fill = "species", facet_grid = strait ~
monsoon) + geom_bar(aes(color = species, fill = species), stat = "identity",
position = "stack") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5,
hjust = 1)) + ggtitle(str_c(one_class, " - Abundant OTUs"))
print(p)
}2.2.9 Heatmaps (abundant OTUs)
for (i in 1:2) {
p <- plot_heatmap(ps_list_abundant$ps[[i]], method = "NMDS", distance = "bray",
taxa.label = "species", taxa.order = "division", sample.label = "sample_label",
sample.order = "sample_label", low = "beige", high = "red", na.value = "beige",
title = ps_list_abundant$title[[i]])
print(p)
}for (i in c(3)) {
p <- plot_heatmap(ps_list_abundant$ps[[i]], method = "NMDS", distance = "bray",
taxa.label = "family", taxa.order = "division", sample.label = "sample_label",
sample.order = "sample_label", low = "beige", high = "red", na.value = "beige",
title = ps_list_abundant$title[[i]])
print(p)
}2.2.10 NMDS
All OTUs
Sample removed because they were pulling the NMDS * PR2X16XS21 it has a single eukaryote (diatom bloom ) * RM13XS36 cause problem for bacteria * PR11XS25 cause problem for hetero euks * SBW11XS26 cause problem for hetero euks * SBW13XS37 cause problem for hetero euks * RM13XS36 cause problem for hetero euks
for (i in 1:3) {
ps_nmds <- ps_list$ps[[i]]
# Remove samples with no reads
ps_nmds <- prune_samples(sample_sums(ps_nmds) > 0, ps_nmds)
# Remove samples that caused problems (1= photo, 2=hetero, 3=prok)
if (i == 1) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR2X16SXS21")),
ps_nmds)
}
if (i == 3) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("RM13XS36")),
ps_nmds)
}
if (i == 2) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR11XS25",
"SBW11XS26", "SBW13XS37", "RM13XS36")), ps_nmds)
}
singa.ord <- ordinate(ps_nmds, "NMDS", "bray")
print(singa.ord)
p <- plot_ordination(ps_nmds, singa.ord, type = "samples", color = "monsoon",
shape = "strait", title = ps_list$title[[i]]) + geom_point(size = 5) +
scale_colour_manual(values = c("red", "orange", "green", "blue")) +
theme_bw()
print(p)
}Square root transformation
Wisconsin double standardization
Run 0 stress 0.17
Run 1 stress 0.18
Run 2 stress 0.17
... New best solution
... Procrustes: rmse 0.058 max resid 0.27
Run 3 stress 0.18
Run 4 stress 0.18
Run 5 stress 0.17
... New best solution
... Procrustes: rmse 0.044 max resid 0.21
Run 6 stress 0.17
... New best solution
... Procrustes: rmse 0.054 max resid 0.21
Run 7 stress 0.18
Run 8 stress 0.18
Run 9 stress 0.17
... New best solution
... Procrustes: rmse 0.047 max resid 0.2
Run 10 stress 0.17
Run 11 stress 0.18
Run 12 stress 0.17
Run 13 stress 0.18
Run 14 stress 0.17
Run 15 stress 0.18
Run 16 stress 0.17
... Procrustes: rmse 0.046 max resid 0.21
Run 17 stress 0.18
Run 18 stress 0.18
Run 19 stress 0.18
Run 20 stress 0.18
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.17
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.12
Run 1 stress 0.12
... New best solution
... Procrustes: rmse 0.091 max resid 0.33
Run 2 stress 0.12
... New best solution
... Procrustes: rmse 0.08 max resid 0.28
Run 3 stress 0.12
... Procrustes: rmse 0.073 max resid 0.22
Run 4 stress 0.12
... New best solution
... Procrustes: rmse 0.084 max resid 0.21
Run 5 stress 0.12
Run 6 stress 0.12
Run 7 stress 0.12
Run 8 stress 0.12
Run 9 stress 0.12
Run 10 stress 0.12
Run 11 stress 0.12
Run 12 stress 0.12
Run 13 stress 0.12
Run 14 stress 0.12
Run 15 stress 0.12
Run 16 stress 0.12
Run 17 stress 0.12
Run 18 stress 0.12
Run 19 stress 0.12
Run 20 stress 0.12
*** No convergence -- monoMDS stopping criteria:
11: no. of iterations >= maxit
9: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.12
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.12
Run 1 stress 0.13
Run 2 stress 0.14
Run 3 stress 0.14
Run 4 stress 0.13
Run 5 stress 0.12
... New best solution
... Procrustes: rmse 0.059 max resid 0.25
Run 6 stress 0.13
Run 7 stress 0.14
Run 8 stress 0.14
Run 9 stress 0.13
Run 10 stress 0.13
Run 11 stress 0.12
... New best solution
... Procrustes: rmse 0.038 max resid 0.26
Run 12 stress 0.12
... Procrustes: rmse 0.038 max resid 0.26
Run 13 stress 0.14
Run 14 stress 0.14
Run 15 stress 0.13
Run 16 stress 0.14
Run 17 stress 0.14
Run 18 stress 0.14
Run 19 stress 0.13
Run 20 stress 0.14
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.12
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Abundant OTUs
for (i in 1:3) {
ps_nmds <- ps_list_abundant$ps[[i]]
# Remove samples with no reads
ps_nmds <- prune_samples(sample_sums(ps_nmds) > 0, ps_nmds)
# Remove samples that caused problems (1= photo, 2=hetero, 3=prok)
if (i == 1) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR2X16SXS21")),
ps_nmds)
}
if (i == 3) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("RM13XS36")),
ps_nmds)
}
if (i == 2) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR11XS25",
"SBW11XS26", "SBW13XS37", "RM13XS36")), ps_nmds)
}
singa.ord <- ordinate(ps_nmds, "NMDS", "bray")
print(singa.ord)
p <- plot_ordination(ps_nmds, singa.ord, type = "samples", color = "monsoon",
shape = "strait", title = ps_list_abundant$title[[i]]) + geom_point(size = 3) +
scale_colour_manual(values = c("red", "orange", "green", "blue")) +
theme_bw()
print(p)
p <- plot_ordination(ps_nmds, singa.ord, type = "taxa", color = "division",
shape = "division", title = ps_list_abundant$title[[i]]) + geom_point(size = 3) +
theme_bw()
print(p)
}Square root transformation
Wisconsin double standardization
Run 0 stress 0.16
Run 1 stress 0.16
... New best solution
... Procrustes: rmse 0.047 max resid 0.23
Run 2 stress 0.17
Run 3 stress 0.17
Run 4 stress 0.16
... New best solution
... Procrustes: rmse 0.041 max resid 0.18
Run 5 stress 0.17
Run 6 stress 0.17
Run 7 stress 0.17
Run 8 stress 0.16
Run 9 stress 0.17
Run 10 stress 0.16
Run 11 stress 0.17
Run 12 stress 0.16
Run 13 stress 0.16
Run 14 stress 0.17
Run 15 stress 0.16
Run 16 stress 0.16
Run 17 stress 0.16
Run 18 stress 0.17
Run 19 stress 0.17
Run 20 stress 0.17
*** No convergence -- monoMDS stopping criteria:
20: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.16
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.093
Run 1 stress 0.092
... New best solution
... Procrustes: rmse 0.041 max resid 0.15
Run 2 stress 0.092
... Procrustes: rmse 0.032 max resid 0.12
Run 3 stress 0.092
... New best solution
... Procrustes: rmse 0.034 max resid 0.13
Run 4 stress 0.092
... New best solution
... Procrustes: rmse 0.023 max resid 0.15
Run 5 stress 0.092
... Procrustes: rmse 0.041 max resid 0.13
Run 6 stress 0.093
Run 7 stress 0.092
... Procrustes: rmse 0.028 max resid 0.13
Run 8 stress 0.092
... New best solution
... Procrustes: rmse 0.038 max resid 0.12
Run 9 stress 0.092
Run 10 stress 0.093
Run 11 stress 0.093
Run 12 stress 0.093
Run 13 stress 0.092
... Procrustes: rmse 0.038 max resid 0.17
Run 14 stress 0.092
... New best solution
... Procrustes: rmse 0.032 max resid 0.13
Run 15 stress 0.093
Run 16 stress 0.093
Run 17 stress 0.093
Run 18 stress 0.092
... Procrustes: rmse 0.038 max resid 0.13
Run 19 stress 0.093
Run 20 stress 0.092
... Procrustes: rmse 0.028 max resid 0.13
*** No convergence -- monoMDS stopping criteria:
8: no. of iterations >= maxit
12: stress ratio > sratmax
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.092
Stress type 1, weak ties
No convergent solutions - best solution after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
Square root transformation
Wisconsin double standardization
Run 0 stress 0.12
Run 1 stress 0.12
Run 2 stress 0.15
Run 3 stress 0.12
Run 4 stress 0.15
Run 5 stress 0.15
Run 6 stress 0.16
Run 7 stress 0.14
Run 8 stress 0.13
Run 9 stress 0.13
Run 10 stress 0.16
Run 11 stress 0.12
Run 12 stress 0.15
Run 13 stress 0.14
Run 14 stress 0.12
... Procrustes: rmse 0.00014 max resid 0.001
... Similar to previous best
Run 15 stress 0.12
Run 16 stress 0.16
Run 17 stress 0.13
Run 18 stress 0.14
Run 19 stress 0.14
Run 20 stress 0.13
*** Solution reached
Call:
metaMDS(comm = veganifyOTU(physeq), distance = distance)
global Multidimensional Scaling using monoMDS
Data: wisconsin(sqrt(veganifyOTU(physeq)))
Distance: bray
Dimensions: 2
Stress: 0.12
Stress type 1, weak ties
Two convergent solutions found after 20 tries
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on 'wisconsin(sqrt(veganifyOTU(physeq)))'
2.2.11 Network analysis
for (i in 1:3) {
ps_nmds <- ps_list_abundant$ps[[i]]
# Remove samples with no reads
ps_nmds <- prune_samples(sample_sums(ps_nmds) > 0, ps_nmds)
# Remove samples that caused problems (1= photo, 2=hetero, 3=prok)
if (i == 1) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR2X16SXS21")),
ps_nmds)
}
if (i == 3) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("RM13XS36")),
ps_nmds)
}
if (i == 2) {
ps_nmds <- prune_samples(!(sample_names(ps_nmds) %in% c("PR11XS25",
"SBW11XS26", "SBW13XS37", "RM13XS36")), ps_nmds)
}
if (i < 3) {
p <- plot_net(ps_nmds, distance = "(A+B-2*J)/(A+B)", type = "taxa",
maxdist = 0.4, color = "class", point_label = "genus") + ggtitle(ps_list_abundant$title[[i]])
} else {
p <- plot_net(ps_nmds, distance = "(A+B-2*J)/(A+B)", type = "taxa",
maxdist = 0.4, color = "class", point_label = "family") + ggtitle(ps_list_abundant$title[[i]])
}
print(p)
}